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Abstract. Wc analyze the problem of the analytical characterization of the 
probability distribution of financial returns in the exponential Ornstein-Uhlenbeck 
model with stochastic volatility. In this model the prices are driven by a Geometric 
Brownian motion, whose diffusion coefficient is expressed through an exponential 
function of an hidden variable Y governed by a mean-reverting process. We derive 
closed-form expressions for the probability distribution and its characteristic function 
in two limit cases. In the first one the fluctuations of Y are larger than the volatility 
normal level, while the second one corresponds to the assumption of a small stationary 
value for the variance of Y. 

Theoretical results are tested numerically by intensive use of Monte Carlo 
simulations. The effectiveness of the analytical predictions is checked via a careful 
analysis of the parameters involved in the numerical implementation of the Euler- 
Maruyama scheme and is tested on a data set of flnancial indexes. In particular, we 
discuss results for the German DAX30 and Dow Jones Euro Stoxx 50, finding a good 
agreement between the empirical data and the theoretical description. 
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1. Introduction and motivation 

Since the pioneering work of Bachelier [1] about tlie fair price of derivatives contracts 
exchanged on the Paris market in the early '900, diffusion processes have been natural 
candidates for the modehng of the stochastic evolution of financial quantities. Few 
years later Einstein and, independently, Smoluchowski published their famous works 
about the diffusion of Brownian particles inside a suspension [21 [3] . Their works were 
later revisited by Langevin [1], but from a completely different point of view. Indeed 
he took the point of view of the single particles and derived a description in terms of 
differential equations governing the microscopic dynamics under the effects of random 
collisions. What is more important for our work, he also focused on the equation 
driving the velocity of the diffusing particles and showed the existence of a stationary 
regime in which the velocity reaches a constant value. He provided the first example 
of what we currently call in a modern language a mean-reverting process described 
by a stochastic differential equation (SDE). Similar arguments but in the language of 
Fokker-Planck partial differential equation can be found in [5l [6] . In financial context it 
is quite common to meet quantities that by construction can not indefinitely diffuse but 
reasonably fluctuate around a stationary value. For example, the modeling of interest 
rate dynamics developed during the '90s entirely deals with mean-reverting rates [8]. 

Empirical studies have also shown that the volatility, the variable that governs 
the amplitude of returns fluctuations, is not constant, as postulated by Black and 
Scholes in their seminal work about option pricing [9] . Indeed, once a suitable volatility 
proxy is defined, it has been empirically demonstrated that it fluctuates along time 
switching between regimes with higher and lower activity (bursting effect). Moreover, 
the assumption of a volatility with a stochastic nature represents a quite effective 
mechanism responsible for the excess of kurtosis observed in the empirical probability 
returns distributions. A mean-reverting dynamics is therefore a natural candidate for the 
modeling of volatility. Indeed in financial literature this idea has been widely exploited, 
as the introduction of several stochastic volatility (SV) models clearly demonstrates 
[ini [HI [121 [13 [in [15] • The research in this field has also been strongly enhanced by 
the pioneering work of Carr and Madan [16] and, more recently, Lewis [17], through 
the introduction of Fast Fourier Transform (FFT) numerical techniques. Also the 
physicists' community has independently analyzed and developed models to capture 
the stochastic nature of volatility, with particular emphasis on the comparison of the 
models predictions with real market data. For a review the reader can see [18], while 
specific analysis are addressed in[l9l[20l[2ll[2a[23l[2a[25l[26l[27l[28]. Our work 
follows in particular the guidelines of [20] and [25]. Both these articles present an 
analytical characterization of the probability distribution of financial returns under the 
assumption of a Geometric Brownian motion coupled through the diffusion coefficient 
with a SDE that describes the evolution of the volatility. Dragillescu - Yakovenko 
[20] and Masoliver - Perello [25] considered two different models, the former authors 
presenting an analysis for the Heston one [13], while the latter focusing on the Scott 
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model [TU]. Main results have been derived under the hypothesis that the stochastic 
process that describes the volatility dynamics has reached a stationary regime, a 
reasonable assumption if the relaxation time is negligible with respect to the time horizon 
along with the process evolves. Following Masoliver and Perello we consider the Scott or 
exponential Ornstein-Uhlenbeck model, in particular for its well known ability to capture 
some stylized facts observed in real market data such as squared-returns autocorrelation, 
leverage effects and multiple time scale properties [22l [25| [26] . However, we extend their 
analysis relaxing the request of the complete volatility thermalization. A good reason for 
our choice is that in many financial applications the involved time horizon is comparable 
with the relaxation time of the process. Moreover, some exotic financial instruments 
do depend on the entire history of the process. The characterization of the probability 
distribution, even if in approximated form, in the out-of-stationary regime can provide 
some interesting insight in the understanding of the process evolution and be of practical 
interest in the field of quantitative finance. 

The structure of the article is the following. In Section [2] we review the SV model we 
will deal with providing the reader with the formulation in terms of systems of stochastic 
differential equations. Section [3] and H] share the same structure. Each of them consists of 
two subsections, the first one depicts the framework and describes the analytical results, 
while the second puts the theoretical results on numerical tests. Firstly we present a 
closed-form expression for the returns probability distribution in the limit case of log- 
volatility fluctuations higher than the volatility normal level. After providing evidences 
that the approximate solution is effective only in the limit of small log- volatility variance, 
we derive an exact closed-formula for the characteristic function but for a linear version 
of the S V model. Again we check the goodness of the results by means of FFT algorithms 
and we test them against the Monte Carlo (MC) simulation of the linear and complete 
models. Section [5] is devoted to the comparison of model predictions with real market 
data. We discuss in detail the calibration procedure and we apply it on a set of financial 
indexes from the equity sector. We report more relevant results for the German DAX30 
index and the Dow Jones Euro Stoxx 50, a market capitalization-weighted index of 50 
European blue-chip stocks from those countries participating in the European Monetary 
Union. The final Section [6] draws the relevant conclusions and suggests some possible 
applications of the analytical results in the field of financial option pricing. Some of the 
results presented in Section [3] are partially shared with a paper [29] appeared on the 
on-line archive during the completion of this work. 

2. The Model 

The process we are investigating is a slightly modified version of the stochastic volatility 
model originally introduced by Scott [10]. In financial and econophysics literature it is 
commonly referred to as an exponential Ornstein-Uhlenbeck process. However, since 
we allow for a non-null stationary mean value for the random variable driving the log- 
volatility of the returns, in our case it should be more correct to speak of an exponential 
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mean-reverting process. A major difference with Scott model is the existence of a 
parameter p, that takes values in [—1, 1] and represents the correlation between the two 
sources of noise in the diffusive process. As we will see later, p different from zero is 
crucial to obtain a non trivial behavior for the skewness of the returns distribution. 
Previous considerations translate in the following coupled SDEs 

dS{t) = fiS{t)dt + me^(*)^(t) dWi{t) , 

S{to) = So , 

dY{t) = a(7 - r(t))dt + kpdWiit) + k^l-pMW2{t) , 
Yito) =yo, 

where dWi and dW2 are two independent Wiener processes, while sq, p, m, yo, a, 7, 
k and p are constant parameters that can be calibrated on real data. If we define 
a{t) = me^^*-* and we set m > in Eq. ([T]) we recognize the same dynamics of a 
Geometric Brownian motion, with a constant drift coefficient p while the so-called 
volatility, characterizes the amplitude of the fluctuations of S. However, the model does 
not only allow for a time varying a, but through the exponential function promotes the 
volatility itself to a random variable driven by the dynamics of Y. By acting on the 
value of > we can strongly modify the behavior of Y. The case k = switches off 
the stochastic nature of Y and it evolves deterministically towards its stationary value 
7, with a characteristic time 1/a (a > 0). If /c is strictly greater than zero, it is well 
known that Y follows a Gaussian process, with 

E [Y] = {yo - 7)e-"(*-*«) + 7 7 (3) 
E [Y^] - E [Yf = — [e-2-(*-*o) + 1] *^ (3 = — (4) 

We have used the notation E [■] to indicate the expectation with respect to the 
probability distribution of Y, so E [Y] and E [Y"^] — E [Y]^ correspond to the usual 
mean value and variance, respectively. For later convenience, Eq. ([1]) can be simplified 
into the following 

dX{t) = -^m^e^^Wdt + me^(*) dWi{t) , (5) 

with boundary condition X{to) = 0, by applying Ito's Lemma to the centered 
logarithmic return X{t) = \n S(t) — \n So — pit — to). 

Equations ([2]) and (I5j) summarize the model we want to focus on. In particular, 
once a time horizon t > to has been fixed, we are interested in the characterization of 
the returns transition probability distribution px{X(t) \Xo, Yo). The returns distribution 
Px{X{t)) can be readily obtained taking the expectation over the initial distribution of 
Xq and Yq. In |25j the authors derive an approximate expression for px under the 
assumption of a stationary regime for the Y variable, i.e. normally distributed with 
mean 7 = and variance /9 > 0, pyq ~ -^(7 = 0,P). The further crucial hypothesis 
assumed by the authors to obtain their results was A = k/m ^ 1. This fact was 
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supported by empirical evidences based on the analysis of daily Dow Jones Index returns. 
In the following section we generalize Masoliver and Perello's result consistently with 
the boundary condition of Eq. ([2]), Fq ~ ^(^o — yo), that is relaxing the request of initial 
stationarity, and allowing for 7 different from zero. We do not relax the condition for 
A in order to follow the same solving strategy of the forward Fokker-Planck equation. 
However we will show the effectiveness of this assumption at the end of the next section 
providing the results of the numerical MC tests. 

Before deriving our results, it can be worth remembering that the distribution of 
returns is surely not the only interesting statistical feature of mean-reverting SV models 
and literature deals with other properties, such as squared-returns autocorrelation, 
leverage effects {i.e. past returns and futures volatilities correlation) and multiple time 
scale properties [25l [26l [30], [SI [32]. However our interest is guided by a possible 
applications of a closed-form expression for the returns distribution in the field of 
financial option pricing. Indeed, in a forthcoming work in preparation [33] , we will show 
how the model under consideration here can emerge in a natural way in a risk-neutral 
description of returns dynamics. 



3. Limit case I: /c/m ^ 1 

3.1. Approximate closed-form expression for the returns distribution 

Given the dynamics ([2]) and ([5]), the associated transition probability density function 
p{x,y\xo,yo) satisfies the forward Fokker-Planck equation [7|,[3l] 

dp ^ m^^,^dp ^ m'^,^d^p ^^dp ^ Jjyp) ^ 1 ^^,d'p ^ ^^^u d^eVp) 
dt 2 dx 2 dx"^ dy dy 2 dy"^ dxdy 
with initial condition 

p{x,y\xQ = 0, yo) = p{x, y\yo) = 5{x)5{y - yo) . (7) 

Instead of working with x, y and t, it is more convenient to introduce the dimensionless 
variables 

T = k'^{t — to), u = \x and v = Xy . 

With respect to r, u and v Eq. ([6]) becomes 

dp _ e^^/^ dp ^ e2^'/^ d'^p a-fXdp^ad (vp) ^X"^ d^p ^ ^d"^ (e^'/^p) 
dr 2X du 2 du^ k"^ dv k"^ dv 2 dv"^ dudv 
with initial condition 

piu,v\yo) = 5{u)5{v - Xyo) . (9) 

The previous relations can be rewritten in terms of the characteristic function 

ip{LU^,U2,T\yo) = f dne*-^" / dv e^'^'M^, v\yo) • (10) 



The condition r = becomes 

V9(c^i,o^2,0|i/o) = / due'^^^6{u) I dve'^'^6{v - Xyo) = e^'^^^^" , (11) 
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while Eq. reduces to 



^(^i,W2,r|?/o) = -^V'(^i,^2-y,r|yo)-^^iV'(^i,^2-y,r|yo) 
H — 'fi^i, ^2, T\yo)-——^—{uJi, UJ2, T\yo) 



fc2 



k'^ duj2 



^2 i 

-Y^2V5(^i'^2,^|yo)-pAu;iCJ2</3(c^i,cJ2--^,7-|yo) • (12) 
In order to solve Eq. ( fT2i) . following [25], we try the ansatz 

Lp{uJi,UJ2,T\yo) = Q-M'^l,yo,r)u:l-B{wuyo,r)w2-C{wuyo,r)+O{w^2) . (13) 

In Appendix A we show how to obtain the ordinary differential equations (ODEs) 
satisfied by A, B and C in the limit of small 002 and up to order 1/A. Once 
the approximate expression of C{LJi,yQ,T) is derived, we can compute the marginal 
distribution pj!(:(x I i/o) noticing that 



Px[x\yo) = — 
Zn 

1 

~ 27r 



-lUllX 



00 
+00 



(p{uJi/X,uJ2 = 0,T\yo)dLJi 



-iuiix^-C[uji/\,yo,T)^^^ 



(14) 



It is worth observing that C corresponds to the logarithm of the characteristic function 
changed by sign. Hence we can easily obtain the cumulant of order n of px 

d^CMX,yo,T 



kn = -{-^y 



(15) 



a;i=0 



This suggests that, instead of computing previous integral directly, we can derive an 
approximated but closed- form solution by exploiting the Edgeworth expansion [35] . We 
report below the expansion considering correction to the Gaussian distribution only up 
to fourth normalized Hermite polynomial [36] 



Px[.x\yo) 



\ (x-kl) 

e 



V2^ 



2fc2 X 



X — kl 



2Akl 



X — kl 



k2 



The explicit expressions of the first four cumulants are given by 

2 

A;2 = — [(l + 27)C + 2(yo-7)(l-e-^)] , 
a 

k3 = [C(l + 7) + {yo - (1 + 27))(1 - 



m'^k'^ 



(16) 

(17) 
(18) 
(19) 



a" 



2C + (1 -e-2^) -4(1 -e" 



+ 4p2 (^ + ^e^f_2(l-e^0) 
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- Vyo (Ce-^ - (1 - e-0 + ^C'e-^) 

+ 4A(^C + 2Ce-^-3(l-e-'^) + ^CV<^ , (20) 

where we have introduced the ancillary variable ( = a(t — to). We have checked the 
consistency of our results when t approaches to with the initial time condition ([7]). 
Indeed, at lowest order in (, ki and /c2 scale linearly, while being k-^ = O(C^) and 
/c4 = 0{(^) the skewness ^ = k^/k^^'^ and kurtosis k = k^^/k^ do not diverge, as expected. 

The analytical expression for C given by Eq. flA.lOp (and from it, by simple 
derivation, the explicit characterization of the cumulants) represents one of the main 
results of this work. However, Eq. fll6p has to be considered with care for various 
reasons. First, it is well known that for an arbitrary choice of /ci, /c2, k^ and k^ the 
positive definiteness of the truncated Edgeworth expansion is not guaranteed. Secondly, 
when dealing with the problem of probability distribution characterization for the sum 
of n random variables {e.g. identically distributed and with finite moment of every 
order) higher order terms of the Edgeworth expansion scales with the power of n"*"/^ 
for suitable r. For example the term proportional to iJa scales with n"^/^, while the 
terms proportional to H4, and (not considered here) scale with n~^. But for the 
case under consideration we have n = 1, so in general we can not a priori estimate the 
goodness of the approximation. Moreover, the entire procedure previously described for 
the derivation of Eq. ( |T6l) is strongly based on some approximations. 

It is worth noting that similar approximations have been adopted in |29j, where 
the complete characteristic function is given by Eqs. (35) and (36) yielding, in their 
notation, the following formula 

f uj"^ 
(p{uj/X, at) = exp< —iLUfi{t) — [rfi^t + 2d{t, 20)] — + ip<i{t, Zo)ur' 

+ {K{t) + d{t, zo) V2) + 0(1/A^) I , (21) 

where the various arguments entering the exponential are defined in [29] . Following the 
procedure of retaining the term inside the exponential function, while performing 
a Taylor expansion starting from the d dependent contribution, the authors of [29] get 
the approximate expression (Eq. (41) in the same paper) 

f i2 1 

-2+'^ I fl ^ ^, ,2 I ■^(+ „ \, ,3 



t—\[l-d{t,Zo 



(p{uj/X, at) = exp< —iujjj,{t) — m t— > 1 1 — i^{t, Zo)uj + ip<^(t, Zo)uj 



+ («:(t) + ^(t, Zo)V2) uj' + 0(1/A^)] . (22) 

This leads to a distribution of returns containing a H2 Hermite polynomial contribution 
as given by Eq. (42) in [29]. 

On the other hand, if one follows a different approach expanding the exponential 
starting from the a;^ term, the approximate characteristic function reads, according to 
the notation of [29], as follows 



^2 ^ 

[rfiH + 2§{t, zo)] Y I [1 + 



Lp(uj/X, at) = exp< —iujjj{t) — [m^t + 2'd{t, Zq)] n 1 + ^P?(^) ^o)^ 



,3 
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+ (Kit) + ^{t, zoY/2) uj^ + 0(1/A5)] . (23) 

This leads to the probabihty distribution of the form of Eq. (fT6l) . where a H2 polynomial 
term is absent by construction and incidentally this expansion corresponds to Edgeworth 
series truncated to the fourth Hermite polynomial [351 EH] . 

In the next section we will discuss numerical tests based on intensive MC simulations 
in order to check the robustness of our results. 

3.2. Numerical Results 

Standard numerical techniques to simulate random paths from the dynamics Eq. ([2]) 
and Eq. ([5]) are widely discussed in literature and a classical reference is given by [37] . 
In particular Chapter 6 is dedicated to discrete schemes of SDEs suitable to generate 
paths using Monte Carlo methods. In our analysis we implement the Euler-Maruyama 
scheme. The two crucial parameters to be chosen correspond to the discrete time step 
At = (t - to)/NSTEP and the total number of MC paths to be generated, MCPATHS. 
The time variables are measured in yearly units, so t — to = 1 corresponds to a diffusion 
process that evolves for one year. We started with At = 10"^ and then decreased it 
until At = 10"^. At the same time we increased MCPATHS from 10^ to 5 ■ 10^. We 
finally fixed At = 10"'^ and MCPATHS = 5 ■ 10^, because we did not find a significant 
improvement with smaller value of At. 

All the numerical MC results have been obtained with a private Micro Beowulf 
Cluster made of 4 nodes, each of which is an AMD Athlon 64 Dual Core with 2.00 GHz 
CPUs and 2.00 GByte of RAM, developed by one of the authors {G.B.) [38]. The resort 
to 64 bit architecture can be useful to avoid some problem that can arise with high 
Monte Carlo statistics. 

In Fig. [1] we report the returns distributions obtained via MC simulation and the 
analytical predictions given by the Edgeworth expansion [161 Curves from bottom to 
top correspond to increasing values of [3: 0.5%, 1%, 2%, 5%, 10%, 25% and 50%, 
p = —0.9, m = 0.1, a = 10, 7 = 0, I/O = 0. From the relation = 2 ■ 10^/5, we derive 
A = 3.16,4.47,6.32,10,14.14,22.36,31.62, that are values consistent with the assumed 
approximation A ^ 1. The choice 7 = guarantees to preserve the interpretation of m 
as the normal level of the volatility in stationary regime. Indeed, 7 7^ introduces an 
exponential correction which we could easily deal with by a suitable redefinition of m 
and a linear shift of Y. The value m = 0.1 implies a yearly volatility of the X process 
of order 10%, while a equal to 10 corresponds to a relaxation time of 0.1 years. In 
Fig. [H we consider the time interval t — to = 0.1 in order to test the goodness of the 
analytical approximation in a non-stationary regime. The agreement is quite good for 
low values of f3 (and therefore low A, but always greater than one) and rapidly worsen 
for higher values for non- stationary regimes. The peak of the MC distribution moves 
rightward and the left tail becomes fatter with respect to the analytical prediction as 
can be clearly seen in the bottom panel of Fig. [H Indeed for the highest f] values the 
Edgeworth expansion starts to oscillate on the tails and becomes negative. In Fig. [2lwe 
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Figure 1. Returns distributions from MC simulation of the Euler-Maruyama scheme 
and analytical probability distributions given by Eq. (|16|) . Parameters values as 
discussed in the text, curves from bottom to top correspond to increasing values of 
/3: 0.5%, 1%, 2%, 5%, 10%, 25% and 50%. Curves have been shifted upwards for the 
sake of readability. Top panel: p = —0.9 and t — to = 0.1; bottom panel: left tails in 
semi-logarithmic scale. 
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show the behavior of the distribution for a positive value of the correlation coefficient, 
p = 0.5, and in stationary regime t — t^ = 1. After ten relaxation times the Y process 
has completely thermalized and again we checked the agreement between MC based 
results and the Edgeworth expansion. As correctly predicted by Eq. (fT9l) . p governs the 
sign of the skewness and the figures confirm that the asymmetry of the distributions is 
opposite with respect to the previous case. Again the best agreement corresponds to 
the lowest /5, while highest values for A are not sufficient to improve it. As in Fig. [H 
significant differences are present in the distribution tails. 

In order to make our considerations more quantitative, we present in Fig. [3] the 
scaling of mean, variance, skewness and kurtosis with time. In each panel we report the 
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theoretical prediction and the numerical MC results with 95% error bars for j3 = 0.5% 
and l3 = 5%. Tab. [1] contains similar results but for all the /? considered in Fig. [T]and 
t-to = l. 

For f3 < 10% mean and variance computed with Eq. (ITTl) and Eq. (fT8|) are 
in statistical agreement with the MC estimates and the linear scaling with time 
theoretically predicted is confirmed by the numerical simulations. The compatibility 
of ^ is limited to small values of j3 values while the situation is quite unsatisfactory 
for K, and (3 > 2%. All these empirical evidences strongly suggest that the expansion 
of Eq. (fT2|) for uj2 nearly zero and up to order 1/A predicts the correct behavior of 
C{uJi,yo,T) for small P values only. The condition A ^ 1 is not enough to guarantee 
that the Edgeworth expansion f|T6l) is a good approximation of the true distribution. 
When the stationary variance of the log-volatility becomes high, the fluctuations of Y 
sensibly deviate from 7 and the exponential function enhances the fluctuations of X. 
This mechanism is responsible for the growth of the empirical kurtosis at high (3. 
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Table 1. Scaling of normalized cumulants for increasing values of /3. The time horizon 
is given by t — io = 1, while the other parameters are the same of Fig. [TJ The index 
refers to values numerically computed with MC simulation (between parenthesis 
we report the error on the last significant digit, 95% confidence level), while cumulants 
with index correspond to Eq. (HH-Eq. ([20]) . 
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0.26 


0.51 


1.29 


2.6 
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4. Limit case II: expansion for small (3 

The numerical results shown in Sect ion [3] suggest that the level of the stationary variance 
of y is a crucial parameter to describe accurately the returns distribution and for this 
reason we focus on it. The intuition behind our choice is that, since (3 governs the level 
of the fluctuations of Y around its stationary value, keeping (3 low allows us to linearize 
the exponential form of the volatility. Moreover we will show in the next subsection how 
it is possible to exactly solve the Fokker-Planck equation associated with the system 
of SDEs by means of linearization. We limit our investigation to < /? < 0.1. Higher 
values can be explored and the goodness of the approximation can be tested comparing 
the numerical results obtained via MC simulation of the linear and complete dynamics. 
However, it is important to note that in practical applications it is usually required to 
keep the volatility non negative. The probability for the volatility to become negative 
can be easily computed and is given by the following formula 

l+7 + (yo-7)e-"(-°)^ — ^Erfcri±jl , (24) 




v//3(e"2-(*-*o) + 1) J 2 V 

where Erfc is the complementary error function. After few relaxation times 1/a and 
for 7 of order we have that the probability reduces to one half of Erfc(\/2a/k). For 
(3 = 1% we have a probability of order 10~^^ , while for /3 = 10% it increases to 4 x 10~^, 
which is still a very small value. 

4.I. Exact solution for the linear model 

The starting point of our analysis is the linearization of Eq. ([5]). Since when the Y 
process termalizes the random variable fluctuates around 7, it is quite natural to expand 
the exponential around the stationary value. By defining m = me' and introducing the 
random variable Z = F — 7 + 1, Eq. ([5]) and Eq. ([2]) can be rewritten as 

— 2 

Tfl 

dX = -^(2Z - l)dt + mZdWi , ^25) 



X(t 



oj 



dZ = a{l-Z)&t + kpmi{t) + k^l- p^m2{t) , ^^g^ 
Zito) = yo - 7 + 1 , 

To derive the analytical expression of px(X(t) |Xo, Zq) we will follow a strategy analogous 
to the one pioneered by Heston in [13]. For notational convenience we indicate 
px{x\xo, zq) shortly as px- The Fokker-Planck backward equation satisfied by px is 
readily written as 
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Heston technique essentially reduces to the observation that, if we assume as a final 
time condition the expression e*^*^, Eq. fl27j) is precisely the partial differential equation 
governing the evolution of the characteristic function f{(j);xo,Zo) implicitly defined by 

/> + oo 

Px{x\xo, ^o) = ^ y e-^-^Vl^; xo, zo)d<P . (28) 
We try a solution of the form 

j^l^^.^j,^^ ^ Q-A{t-to,4>)+B{t-to,<t>)zo+C{t-to,<l,)z^+i(l,xo ^ ^29) 

We substitute it in Eq. ([27]) and we get 

— 2 — 2 

A + Bzo + Czl = — (2zo - l)i4> - —zl{-4>^) - a{l - zo){B + 2Czo) 

[{B + 2Czof + 2C] - pkmzoi(j){B + 2Czo) , (30) 

2 

where the dot stays for a derivative w.r.t. to- If we collect the quadratic and linear 
terms, independent of zq, we have the ODEs satisfied by C, B and A respectively: 

C = —0^ + 2aC - 2k'^C^ - pkThi(j)2C , (31) 

B = rfiHcp - 2aC + aB - 2k'^BC - pkmupB , (32) 

A=-—t^-aB^-iB' + 2C) , (33) 

with final time conditions 

C(O,0) = 0, (34) 
i3(O,0)=O, (35) 
^(0,0) = 0. (36) 

Equation (13T]) is a Riccati type ODE and once it has been solved, we insert the solution 
in Eq. (1321) and Eq. (1331) and we integrate them out in the usual way. The explicit 
expression for C, B and A are quite involved. To improve the readability we define some 
auxiliary variables 

d = 2\/ k'^w?^'^ + (a — ikrhp(j)y , b = 2{a — ikrhpcp) , 

g = , h = ifh^cj) and n = —r^ib — d) . 

b + d 2k^ 

Now the desired functions read 

h — ri 1 — p-d{t~to) 

e--^<i(^-''){{g+l)h-2n)+n + e-^('-'^\n-gh)-h 
^(^ - = 2 rf(l-,e-'^(^~M) ' ^38) 
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and 



A{t - to, , 



2a 



1 



h{t-to) 
{g + l)h-2n 



ln(l - 2 



+ 



-ln(l + V^e-t(*-*"))+ln(l + v/^) 
n{g + 1) -2gh ^^^^^ „„„d(t_to) 



in(i - Va) 



[ln(l - 5(6' 



ln(l - g 



n — h , 
+ ^^{t-to) 



2^ 



d 

{n - ghf 



-d{t-to) _ I 



d^g L(l-(7)(l-^7e-<^(*-*°)) 9 
-\n{l-g) 



+ -i\n{l-ge 



-d{t-to)^ 



{{g + l)h - 2nf + 2{n - gh){n - h) 



~d{t~to) _ I 



rf3 



(1 -5)(1 -5e-'^(*-*o)) 



+ 



(n - hy 
d^ 



{t -t,) + - (ln(l - ^e-'^^*-*")) - ln(l - g)) 



+ 



d{l - g){l - ge-'^(^-^o)) 
{{g+l)h-2n)^ {{g + l)h - 2n){n - h) 



X 



ln(l + 



d^V9 

ln(l + ln(l - V^e- 



H*-*o)^ 



+ ln(l-v/^)^ 
'{{g + l)h-2nf {{g + l)h - 2n){n - h) 



d'g 



d^ 



lit-to) 



)(e- 



-i(t-to) 



[l-g){l-ge- 



h-d 



{t - to) + 



9 



[HI 



9^ 



-d(t-to)^ 



Hl-g)] 



(39) 



4 " ■ 

The three functions A, B and C require some care. Indeed, similarly to the Heston case, 
non trivial problems emerge due to the multi-valued nature of complex square root 
and logarithm. In particular, due to "branching" effects, the characteristic function 
can become discontinuous. Following the same arguments discussed in [39| BOl HI] , we 
checked the smoothness of xo, zo) for the set of parameters used in the next Section. 

The complete characterization of /(</>; Xq, Zo) is the main result of this paper. The 
only remaining task to be performed to obtain the probability distribution in the direct 
space is to Fourier anti-transform it. We stress again that the obtained solution is an 
exact one but for the linear model. In the following Section we provide a set of (3 values 
for which we have numerically tested the effectiveness of the approximation. 
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4-2. Further Numerical Results 

The exact solution for the probabihty distribution of returns in the framework of the 
hnear problem is given in an implicit way by means of a Fourier anti-transformation of 
Eq. (!28l) . which can not be computed analytically. We compared two different ways for 
the numerical evaluation oi px'- the integration with a trapezoidal approximation and 
the FFT algorithm. In order to improve the efficiency of the methods one can observe 
that Px is a real function, so the integrand of Eq. ( l28i) must have its real part even and 
its imaginary part odd. We can integrate on one half of the entire real axis and take 
twice the real part of the result. To sum up: 



We checked the convergence of the two methods mentioned previously increasing 
the number of points for the sampling of the integrand. All the results presented in the 
next figures have been obtained with the FFT algorithm with 2^^ points on the fixed 
interval of integration [0, 10^]. The tails of the resulting distribution must decrease to 
zero, but they are extremely sensitive to the FFT frequency sampling so they have been 
neglected in Fig. |H For the reasons explained in Section 13.21 we chose At = 10~^ and 
MCPATHS= 5 ■ 10^ to simulate the dynamics of the exponential and linear models. 

In Fig. m we report the histograms of returns obtained via MC simulations of the 
exponential and the linear models with the px curve given by the numerical Fourier 
anti-transformation, m = 0.1, a = 10, 7 = and Vo = are fixed for all the figures 
and table of this Section. Each group of three curves corresponds to (3 values 0.5%, 
1%, 2%, 5% and 10%, from bottom to top respectively. Again each curve has been 
shifted for the sake of readability. The top frame has been obtained with a positive p 
value (0.5), which is evident from the rightward asymmetry of the curves, while, in the 
bottom panel, curves have large leftward asymmetry corresponding to p = —0.9. The 
agreement between theory prediction and MC simulation is very good for f] = 0.5%, 
1% and 2%, in both cases of p values. The distributions of returns from exponential 
model are well reproduced by the ones from linear model and, as a consequence, by 
the semi- analytical approximation px- The agreement worsen when (3 increases to 5% 
and 10%, as it can be seen observing the fatter tail in each panel. All the curves in 
Fig. m show that the returns distribution in the linear model is exactly the one computed 
via FFT, thus confirming the goodness of the exact solution. Finally we observe the 
good behavior of all the semi-analytical curves, which, since they are exact solutions of 
a Fokker-Planck equation, do not become negative or fiuctuate like in the case of the 
Edgeworth expansion considered in the previous Section. 

Again, to make the considerations more quantitative, in Fig. [5] the scaling with 
time of mean, variance, skewness and kurtosis is shown for (3 = 0.5% and 5%, p = —0.9. 
Then in Tab. [21 the numerical values of normalized cumulants computed with MC 
samples of returns from exponential and linear models are shown for (3 = 1%, 2% and 
t — to = 0.01,0.1,0.2,0.5,1. The results show that, for /5 < 2%, mean, skewness 




(40) 
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Figure 4. Returns distributions from MC simulation of exponential and linear 
models in comparison with the probability distributions given by Eq. (|28|) . numerically 
computed by FFT algorithm. Parameters values as discussed in the text, curves from 
bottom to top correspond to increasing values of /?: 0.5%, 1%, 2%, 5% and 10%. 
Curves have been shifted upwards. Top panel: p = 0.5 and t — to ~ 1] bottom panel: 
p = —0.9 and t — to = 1. 
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and kurtosis from the complete and linear models agree within the usual 95% MC 
confidence level, while variance values are not statistically compatible. This is an 
expected result, since the linear model can not exactly reproduce the complete one 
and for this reason we also expect that for higher values of MCPATHS the discrepancies 
for the other cumulants should emerge. However, we are able to evaluate the relative 
differences among cumulants values and we can evaluate the degree of approximation 
of the analytical solution: for (3 = 1% the relative disagreement between k^^^ and 
increases with time from 0.6% to 1%. As expected, the disagreement between cumulants 
computed for the exponential and linear models increases with the value of (3. 

In comparison with the results given in Section [321 one can note that the tails of the 
distribution are better reproduced by the semi-analytical approximation (for mentioned 
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Figure 5. From top left clockwise: scaling with time of mean, variance, kurtosis and 
skewness. Comparison between numerically estimated values of normalized cumulants 
for the exponential (Eq. Eq. (O) and the linear models (with MC 95% confidence 
level) . p = —0.9, (3 = 0.5%, 5% and other parameters as in the text. 
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values of /?) than by the Edgeworth expansion, the other parameters being fixed. Also 
the general agreement of cumulants for (3 < 2% is better than the corresponding results 
given in Fig. [3] and Tab. [H 

5. Real data analysis 

We test the effectiveness of previous analytical results on a data set composed by several 
financial indexes. In particular we consider the following indexes from the equity sector, 
DAX30, CAC40, FTSEIOO, S&PMib, S&P500, DJ Euro Stoxx 50 and NYSE, but we 
detail the analysis only for DAX30 and DJ Euro Stoxx 50. Similar results have been 
found also for the other series. 

DAX30 time series is made of 12173 daily close prices, from 4*^^ January 1960 until 
30*^^ June 2008, while for DJ Euro Stoxx 50 we consider 5566 prices, from 1'^* January 
1987 until 3P* August 2008. Time is measured on a yearly base, so for DAX30 and 
DJ Euro Stoxx 50 we have 48.5 years and 21.75 years time windows respectively. The 
model described by Eqs. ([H) and ([2]) depends on 8 parameters. However, sq corresponds 
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Table 2. Scaling with time of normalized cumulants for (3 = 1% and 2%, p = —0.9 
and other parameters as in the text. The indexes and '^"^ refer to values 
numerically computed with MC simulations with exponential and linear dynamics 
(between parenthesis we report the error on the last significant digit, 95% confidence 
level). 









/3 = 1% 






t - to(yr) 


0.01 


0.1 


0.2 


0.5 


1 


^Lin (10-4) 


-0.50(8) 
-0.50(8) 


-6.0(2) 
-4.9(2) 


-9.9(4) 
-9.8(4) 


-25.3(6) 
-24.9(6) 


-50.8(9) 
-49.8(9) 


fc^"P(10-4) 

/cfr'" (10-4) 


1.002(1) 
1.001(1) 


10.15(1) 
10.09(1) 


20.42(2) 
20.26(2) 


51.26(6) 
50.79(6) 


102.5(1) 
101.5(1) 


^Exp 


-0.107(3) 
-0.107(3) 


-0.282(4) 
-0.279(4) 


-0.312(4) 
-0.306(4) 


-0.276(4) 
-0.271(4) 


-0.219(4) 
-0.215(4) 


^Exp 
^Lin 


0.02(1) 
0.02(2) 


0.15(2) 
0.11(2) 


0.18(2) 
0.14(2) 


0.14(2) 
0.11(2) 


0.08(2) 
0.07(2) 








/? = 2% 






/fcf"P(10-4) 
f.Un (iQ-4) 


-0.50(8) 
-0.50(8) 


-5.0(2) 
-4.9(2) 


-10.1(4) 
-9.8(4) 


-25.8(6) 
-24.9(6) 


-51.8(9) 
-49.8(9) 


^Exp(lQ_4) 
/C^'" (10-4) 


1.004(1) 
1.002(1) 


10.28(1) 
10.16(1) 


20.77(2) 
20.45(2) 


52.35(7) 
51.38(6) 


104.8(1) 
102.8(1) 


^Exp 


-0.151(4) 
-0.150(4) 


-0.402(4) 
-0.392(4) 


-0.443(4) 
-0.429(4) 


-0.393(4) 
-0.380(2) 


-0.311(4) 
-0.300(4) 


^Exp 
^Lin 


0.04(2) 
0.03(2) 


0.30(2) 
0.22(2) 


0.36(2) 
0.27(2) 


0.28(2) 
0.21(2) 


0.17(2) 
0.13(2) 



to the asset spot price so the true free parameters are /i, m, yo, a, 7, l3 and p. In order 
to estimate their values we adopt the following strategies: 

11: from the discretized version of equation ([1]) 

^ = fiAt + mVAte^'et, (41) 

where ASi = Si^i — Si and ej ~ A/'(0, 1), we can conclude that the expectation 
{ASi/ Si) /At over the real data sample provides an estimate of fi. For DAX30 
At = 3.98 X 10~^ while for DJ Euro Stoxx 50 we have At = 3.91 x 10^^ 

7, I/O- Remembering the definition cr(t) = me^^^\ it is readily proved that 

with E[y(t)] and E [F(t)^] given in Eqs. (j3]) and (jlj). It is common practice 
to assume that, for the observed time series, the volatility process has already 
reached the stationary state. Under this assumption, previous expression reduces 
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to E = (me^)^ . It is crucial noticing that all the moments do not depend 

on and e''' is always coupled with m. For this reason we introduce the parameter 
fh = me^ (see also lines before Eq.Q) and we set uq equal to zero. 

m,/3: These two parameters completely specify the stationary log-normal distribution 
of the stochastic variable a. To extract the distribution of the hidden variable a 
from the series of financial returns, we implement the methodology described in 
Appendix B of reference [12] . The consistency of our code has been tested over the 
NYSE time series and we have found results in full agreement with those quoted in 

m 

a,p: Finally, to estimate a and p, we search for values able to reproduce the empirical 
scaling with time of real data skewness and kurtosis. We consider time horizons 
from one day to one hundred days and normalized cumulants are evaluated with 
standard estimators. By means of centered returns, computed from market prices, 
we obtain the empirical skewness, <;ph, and kurtosis, hiph, and corresponding errors, 
e''pf^ and epj^. By generating 10000 paths from the exponential Ornstein-Uhlenbeck 
dynamics, we can compute the MC estimators <;mc, i^mc and associated errors, 
e^(^ and elj^. The optimal a and p are given by those values that minimize the 
sum over 100 time horizons of the normalized squared differences, according to the 
following formula: 



100 



(a*, p*) = min > 

a>0,pe(-l,l)^ 

1=1 



(43) 



The subscript i means that the corresponding quantity is evaluated at time iAt. 
Since the presence of noisy denominators in Eq. (143|1 . we can not resort to 
optimization algorithms based on gradient methods. For this reason we implement 
the principal axis approach with the one dimensional search based on the Brent 
method (see opt/ directory at http://www.netlib.org). 

The probability distribution P(cr) of volatility is plotted in Fig. El both for the DAX30 
index and the DJ Euro Stoxx 50. Distributions tails suffer low statistics effects, but in 
the central region they are well fitted by a log-normal distribution 

I \ 1 1 /^ logg-logffo y 
m) = /— exp-- . (44) 



/27rs(T 2 

The fit is performed in the range 0.0004 < cr < 0.015 for DAX30 and gives loguo = 
-4.492 ±0.001 and s = 0.334 ±0.001, while for the DJ Euro Stoxx 50 the log-normal fit 
is consistent in a wider region 0.0015 < a < 0.015 and gives logao = —4.7049 ± 0.0008 
and s = 0.5147 ± 0.0007. Tab. [3] details the corresponding values for fh = ao/\/'Ai 
and /5 = s^. In Fig. [7] we plot the scaling with time of skewness and kurtosis for the 
DJ Euro Stoxx 50 index. We have found similar results for DAX30, but we do not 
report them here. Error bars represent the 68% confidence level (CL), while the solid 
lines correspond to the analytical expressions given by Eqs. (fT9l) and ( l20l) . These lines 
have been generated using parameters values obtained by the minimization of r.h.s. of 
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Figure 6. Probability distribution P{<t) of volatility for DAX30 (blue boxes) and DJ 
Euro Stoxx 50 (violet diamonds). Solid lines correspond to a log- normal fit. 




0.0025 0.005 0.0075 0.01 0.0125 0.015 



Table 3. Estimated values for the model parameters for DAX30 and DJ Euro Stoxx 
50 indexes. 





/"(yr^^) 


yo 


7 


m (yr 2 ) 


/3 


a(yr-i) 


P 


DAX30 


7.39 X 10-2 








14.52 X 10-2 


11.16 X 10-2 


30.76 


-0.54 


DJ Euro Stoxx 50 


7.97 X 10-2 








14.45 X 10-2 


26.51 X 10-2 


47.74 


-0.52 



Figure 7. Scaling with time of skewness (left panel) and kurtosis (right panel) for DJ 
Euro Stoxx 50 index. We plot the scaling of empirical cumulants (violet crossed-line) 
with the corresponding 68% CL bars and cumulants with 68% CL bars from the MC 
simulation (blue boxed-line). The solid line corresponds to the analytical expressions 
given in Eqs. ^ and (^0)) . 
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Eq. (H3|) . Their values are detailed in Tab. [31 It is worth to comment that in order to 
capture the leftward asymmetry of the real data distribution the correlation coefficient 
p is correctly predicted to be negative. The relaxation time 1/a ~ 5 days implies a 
quite fast thermalization process. The plotted results indicate that the exponential 
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Ornstein-Uhlenbeck model is unable to capture the excess of kurtosis observed in many 
high-frequency returns distributions [HlllS], but it can provide reliable predictions for 
returns distributions corresponding to sufficiently long-time lags and is able to account 
for the transition from a leptokurtic to a Gaussian regime, as observed in other SV 
models [201 El]. Finally, we present the comparison between the probability distributions 



Figure 8. DAX30 index returns distributions for two different time horizons, 25 



trading days (first row) and 45 trading days (second row). 




for DAX30, see Fig. [8], and for DJ Euro Stoxx 50, see Fig. [9], computed according to 
the various models we have discussed in this work. First row of both figures is obtained 
setting the time horizon equal to 25 trading days, while the second one corresponds to 
the 45 trading days horizon. The distributions are presented both in log-linear and linear 
scales, in order to allow the reader to appreciate the behavior on the tails and in the 
central region, respectively. The boxed-lines represent the empirical histograms (all the 
bins contain at least ten points), while black points correspond to the histograms from 
the MC simulation of the exponential Ornstein-Uhlenbeck model. Moreover, we report 
the Fourier transform of the characteristic function for the linear model (solid black 
line), the Normal Maximum Likelihood fit (solid red line) and the analytical solution 
(fT6l) (solid blue line). The Normal approximation is scarcely representative of the true 
empirical distribution, while the exponential model captures in a quite effective way 
the leftward asymmetry, the fatter tails and the narrower central region. The Fourier 
transform of Eq. ( l29i) and Eq. ( |T6l) are both candidates to approximate the distribution 
of the exponential model. For the DAX30 index, the level of (3 equal to 11.16% allow 
us to be confident in a good performance of the linear model. This is indeed the case. 
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as it is clearly shown in Fig. [HI The situation slightly worsen for the D J Euro Stoxx 50 
index, where /? = 26.51%. However, also in this case, the linear model performs better 
than the solution based on the Edgeworth expansion, that presents a narrower central 
peak and thinner tails. Moreover, for the r.h.s. of Eq. (fT6!) positive definiteness is not 
guaranteed. This is confirmed looking at the right tails in Fig. [8], where the distributions 
in log-linear scale have to be truncated since the presence of negative values. 



Figure 9. DJ Euro Stoxx 50 index returns distributions two different time horizons, 
25 trading days (first row) and 45 trading days (second row). 




0.001 



Before concluding, it should be emphasized that further considerations about the 
capability of the linear model to reproduce known financial evidences could be derived 
by comparing the theoretical predictions of the model for the fair prices of options with 
the corresponding market option values. Actually, since our model has been mainly 
developed having in mind the time evolution of risk neutral distributions, this would 
allow to derive implied distribution functions over a given time horizon to be compared 
with the expectations of the model. For example, in [T8l 1^ . it is shown that the 
implied distribution functions corresponding to options with a time to maturity of one 
month have an empirical kurtosis value of the order 1, in good agreement with the 
model predictions for f3 values in the range 5% -r- 10%. This kind of analysis, which 
could be performed along the lines of [15], is presently in progress and is left to a future 
publication. 
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6. Conclusions and perspectives 

In this paper we have deah with the problem of the analytical characterization of the 
probability distribution of returns under an exponential Ornstein-Uhlenbeck dynamics 
for volatility. Such processes have been widely studied, in particular in the econophysics 
literature, from the point of view of multi-time scale properties, leverage effect and 
stationary volatility distribution [HI [201 HD 1221 [251 [21]. However, to the best of our 
knowledge, a systematic approach to the study of returns distribution is missing in 
literature. Our interest is referred to the possibility of a financial application in the 
context of option pricing of a closed-form expression, even if approximated, for the 
probability distribution. The first attempt to obtain such a result can be found in 
[25] under the hypothesis discussed in Section [2l in particular under the assumption 
of a stationary regime for the hidden Y process. We have generalized that result to 
the out-of-stationary regime, with the same constraint of the amplitude of volatility 
fluctuations, k, higher than its normal level, m (A = k/m ^ 1). Indeed, we have 
provided a closed-form expression in terms of the Edgeworth expansion. This is one of 
the main results of our work. We have strongly tested our analytical result with MC 
simulations of the discrete Euler-Maruyama scheme of Eq. ([2j) and Eq. ([5]). The goodness 
of the approximation has been checked with a careful analysis of the discretization step 
and MC number of paths for a quite reasonable choice of model parameters. We have 
found a good agreement for low values of /3, the stationary variance of the log-volatility, 
while our results worsen when (3 increases even if A increases too. For this reason 
in Section [H we have explored the scenario opened by low jS values which allows us 
to expand the exponentials involved in the model up to first order in Y. Following 
the technique pioneered by Heston in [13], we have solved exactly the Fokker-Planck 
backward equation associated with the linear model by means of a suitable trial solution. 
The full expression is reported in Section [Hand is the main result of this work. Eq. ([371) - 
Eq. ( !39l) are quite involved and require some cares, as stressed in the text. Indeed 
they naturally arise in a complex domain and the presence of multi-valued square root 
and logarithmic function introduces pain and angers of branching effects. They have 
been tackled following suggestions coming from the literature related to the Heston 
model. In particular, we have verified the smoothness of the real and imaginary part 
of the characteristic function for the parameters we consider. The final task of Fourier 
anti-transformation has been accomplished by means of FFT techniques. We have 
numerically tested the effectiveness of the analytical expression and we have found a 
perfect agreement. Moreover, the agreement between the results of the MC simulation 
of the discrete scheme of both the complete and linear model make us confident of the 
numerical convergence of the exponential process for low (3 values. The last section of 
this paper has been devoted to the comparison of analytical predictions with real market 
data. The results of our analysis, in particular those for the German DAX30 and the 
Dow Jones Euro Stoxx 50 indexes, confirm the capability of the linear model to capture 
the statistical properties of the returns distribution for moderate values of j3 and over 
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appropriate time horizons. 

As future development we plan to show how the SV model we have studied in 
this work can emerge in a quite natural way in the context of option pricing and risk 
management. In particular, we will discuss how to construct a hedged portfolio with 
associated underlying's dynamics following the exponential Ornstein-Uhlenbeck model. 
The knowledge of a closed-form expression for the characteristic function allows to 
implement a Carr - Madan-like approach and to calibrate the Heston-like option price 
formulae over the implied volatility surfaces. Moreover, we are interested in performing 
a numerical comparison between option prices obtained through full MC evaluation, our 
analytical formulae based on a linear model and other closed-form expressions available 
in literature [291 H6] . 
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Appendix A. Derivation of A{ui,yo,T), B{ui,yo,T) and C(c<Ji, yo, t) 

In order to derive approximate relations for the functions A, B and C, we substitute 
the ansatz ( |T3l) in Eq. ( |T2i) . From relation ( |T4l) we argue that relevant information 
correspond to small UJ2 regime and moreover we look for a solution valid in large A 
regime. Expanding the exponentials up to order 1/A and equating the coefficients of uj2, 
UJ2 and terms independent of 0^2, we find the following first order ODEs 

A = -2e(uii)A+\>f , (A.l) 

C^^B,^.^, (A.3, 

with 

^(^1) = ^ - ip^i ■ (A.4) 

The dot denotes a derivative w.r.t. r. The initial condition ffTTl) traduces in 

AK,yo,0) = 0, (A.5) 

5(u;i,?/o,0) = - iXyo , (A. 6) 

CK,yo,0) = 0. (A.7) 

Solutions of equations (lA.ip . (]A.2p and (1A.3|) with boundary conditions (lA.Sp . (]A.6P 
and (1A.7P read respectively 
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C{uJi,yo,T) = -^^ + Y^ 
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